import netCDF4 
import sys, struct
import math as pm
import numpy as np
import pylab as p
import math
import csv
import matplotlib.pyplot as plt
import h5py
from netcdf_reader import *
from scipy import ndimage
from sklearn.decomposition import PCA
#http://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html

FILE_NAME = 'pe_dif_sep2_98.nc' 
FILE_NAME_CENTRAL_FORECAST = 'pe_fct_aug25_sep2.nc'
INPUT_DATA_DIR = '/home/behollis/DATA/pierre/ocean/'
OUTPUT_DATA_DIR = '/home/behollis/Dropbox/visweek2015/oslines/'
#realizations file 
pe_dif_sep2_98_file = INPUT_DATA_DIR + FILE_NAME
pe_fct_aug25_sep2_file = INPUT_DATA_DIR + FILE_NAME_CENTRAL_FORECAST 

TOTMEMS = 16
X = 152
Y = 152
DIM = 2
 
def writeElapsedTime(start, end):
    import os
    systime = os.times()
    tout = open(OUTPUT_DATA_DIR+'timingFtvaOcean.{0}.txt'.format(systime[4]), 'w') #sys time
    tout.write('Elapsed time: {0}'.format(end - start))
    tout.close() 

if __name__ == '__main__':
    
    #realizations reader 
    rreader = NetcdfReader(pe_dif_sep2_98_file)
    landv = rreader.readVarArray('landv')
    landv = ndimage.rotate(landv, -90)
    
    import time
    start = time.time()
    
    output = h5py.File( OUTPUT_DATA_DIR+'stirringFtvaBackward.hdf5', 'w')
    dcov = output.create_dataset(name='cov', shape=(X,Y, 2, 2), dtype='f')
    #dcmps2 = output.create_dataset(name='cmps2', shape=(X,Y, 2, 2), dtype='f')
    #dmean = output.create_dataset(name='mean', shape=(X,Y,2), dtype='f')
    #dncmps = output.create_dataset(name='ncmps', shape=(X,Y), dtype='f')
    #dexpvarratio = output.create_dataset(name='expvarratio', shape=(X,Y,2), dtype='f')
    
    f = h5py.File(OUTPUT_DATA_DIR + 'industrialStirringSlinesTS050.500steps.ss05.2.hdf5', 'r')
    mems = [ d for d in f ]
    
    for x in range( 0, X ): 
        for y in range( 0, Y ):
            tpts = list()
            print 'calculating ftva for: {0}, {1}'.format(x, y)
    
            for mem in mems:
                dir = mem + '/x' + str(x).zfill(3) + '/y' + str(y).zfill(3)
              
                try:
                    #pt0_0 = f[dir][0][0] #end point in forward integration
                    #pt1_0 = f[dir][1][0]
                    
                    pt0_1 = f[dir][0][-1] #end point in backward integration
                    pt1_1 = f[dir][1][-1]
                    
                    #if math.isnan(pt0_0) is not True and math.isnan(pt1_0) is not True:
                    #    tpts.append([pt0_0, pt1_0])
                    if math.isnan(pt0_1) is not True and math.isnan(pt1_1) is not True:
                        tpts.append([pt0_1, pt1_1])
                    
                except:
                    print 'could not read values from hdf5 file'
                    print '\t for x{0} y{1} {2}'.format( x, y, mem )
                    continue
                
                
            if len(tpts) < 2:
                print 'no pca calculated for land mask,etc. '
                continue  
            
            #get PCA for this step...
            '''    
            pca = PCA(n_components=2)
            pca.fit_transform(tpts)
            '''
            
            dcov[x,y] = np.cov(np.array(tpts).T)#pca.get_covariance()
            '''
            dcmps2[x,y] = pca.components_
            dmean[x,y] = pca.mean_
            dncmps[x,y] = pca.n_components_
            dexpvarratio[x,y] = pca.explained_variance_ratio_
            '''
            
            
    end = time.time()
    writeElapsedTime(start, end)

    f.close()
    output.close()
    
    print 'finished!'